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The study of models with extended Higgs sectors requires to minimize the corresponding Higgs 
potentials, which is in general very difficult. Here, we apply a recently developed method, called 
numerical polynomial homotopy continuation (NPHC), which guarantees to find all the stationary 
points of the Higgs potentials with polynomial-like nonlinearity. The detection of all stationary 
points reveals the structure of the potential with maxima, metastable minima, saddle points besides 
the global minimum. We apply the NPHC method to the most general Higgs potential having two 
complex Higgs-boson doublets and up to five real Higgs-boson singlets. Moreover the method is 
applicable to even more involved potentials. Hence the NPHC method allows to go far beyond the 
t— 1 limits of the Grobner basis approach. 

o 
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I. INTRODUCTION 

The Standard Model (SM) comes with a simple Higgs potential: due to gauge invariance and renormalizability, 
there is only one Higgs doublet which appears in the potential in one bilinear and one quartic term. In particular, 
finding the global minimum of the Higgs potential, which gives the vacuum-expectation value of the Higgs field at the 

Oh stable vacuum, is quite straightforward. 

f^ The situation is very different in models with extended Higgs sectors. For instance in the general two-Higgs-doublet 

model (THDM), introduced by T.D. Lee [T] in order to allow for CP violation in the Higgs sector, we encounter 
already three bilinear terms as well as seven quartic terms in the potential, corresponding to 14 real parameters. 
For some recent works on the THDM we refer to [2HB]. Much motivation for the two-Higgs doublet model arises 
from supersymmetry, since supersymmetric extensions of the SM require to have at least two Higgs doublets - 

£>. besides the corresponding superpartners. The minimal supersymmetric extension of the Standard Model (MSSM) 
Q\ has exactly two Higgs doublets and is therefore a special THDM; for a review see |7]. However the MSSM has a 
number of drawbacks, in particular the so-called /i parameter in the superpotential has to be adjusted by hand to 
the electroweak scale, that is, the electroweak scale appears not via electroweak symmetry breaking. This /i problem 
may be circumvented by a vacuum-expectation value of an additional (complex) Higgs-boson singlet, leading to the 
next-to- minimal supersymmetric standard model (NMSSM); for reviews see (5JI9]. Different supersymmetric models 
have been proposed which have two Higgs-boson doublets and one or several Higgs-boson singlets; an overview can 
<N be found in HDl. 



X 



Other examples of extended Higgs sectors appear in the study of additional discrete symmetries. Extending the 
Higgs sector and assigning all fermions and scalars to certain irreducible representations with respect to the discrete 
symmetry (for instance the quaternion group with eight elements, Qs), the masses and mixing angles of the quarks 
C^ and leptons may be given by very few parameters; see for instance [11U12J . 

While dealing with an extended Higgs sector it is mandatory to minimize the Higgs potential. The determination 
of the global minimum is essential to determine the masses of the fermions via the Yukawa couplings and to ensure 
that the electroweak symmetry breaking behavior follows SU(2)l x U(1)y — > U(l) em . However, finding the global 
minimum of a given non-linear multivariate function is a highly non-trivial task. There are a number of numerical 
methods, e.g., conjugate gradient, steepest descent, simulated annealing, genetic algorithm, etc. which are devised 
to find the global minimum but for one reason or another almost all of them are known to fail even for a relatively 
small system: the multivariate functions usually come with many minima and the above mentioned methods can be 
easily stuck to a local minimum instead of the global minimum. Hence, one can never be sure that the actual global 
minimum is found by these methods. 
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A more systematic approach to minimize the potential is given by the solution of the stationarity equations. Suppose 
the potential is bounded from below, then the stationary points with the lowest potential value are the global minima, 
that is, the vacua. However, again, solving the stationary equations amounts to finding the solutions of a system of 
non-linear multivariate equations which is equally hard. However, the finding of all stationary points not only reveals 
the global minimum, but gives some valuable insight to the structure of the potential: in this way all maxima, saddle 
points, local minima besides the global minima are detected. 

If the potential and hence the system of stationary equations have polynomial-like non-linearity (as opposed to 
involving transcendental functions etc.), then one viable option is to use a symbolic method called the Grobner basis 
approach. The Grobner basis approach was used to find all the stationary points and hence the global minimum of 
extended Higgs potentials in [13]. There, for a special case with two doublets and two real singlets, it was shown that 
extended Higgs potentials may develop a very rich structure of stationary points. 

The Grobner basis approach, however, has several severe drawbacks: firstly, the Buchberger algorithm is known to 
suffer from the exponential space complexity, i.e., the memory (RAM) required by the machine exponentially blows 
up with increasing number of variables, number of equations, degree of equations or number of monomials. The 
Buchberger algorithm is highly sequential and hence very difficult to parallelize. Thus, practically, the Grobner basis 
approach is limited to only small polynomial systems of equations of lower degrees. 

In particular, in the example of two doublets and n real singlets it is currently not possible to go far beyond the 
case n = 2, as considered in |13j . 

In this work we shall introduce the numerical polynomial homotopy continuation method |14l 1 15] as a method to 
find all stationary points of extended Higgs potentials. We shall demonstrate that this method is very powerful and 
overcomes all of the above mentioned shortcomings of the Grobner basis approach. The approach allows to identify 
all stationary solutions and supposed the potential is bounded from below the global minima arc therefore certainly 
detected. As an example we minimize the most general potential with two complex Higgs-boson doublets and n real 
Higgs-boson singlets. We do a detailed analysis for the cases up to n — 5 in this paper, but we emphasis that we can 
straightforwardly go up to n — 10 even on a regular desktop machine. And by exploiting the parallelizable nature of 
the method, it is possible to go even beyond the n = 10 case. 

II. THE HIGGS POTENTIAL AND THE STATIONARITY EQUATIONS 

In this Section, besides fixing the notation we explicitly write down the potential and its stationary equations along 
with the tadpole conditions. 

Let us consider a Higgs potential with two Higgs-boson doublets and n Higgs-boson singlets [T3]. By convention 
we assume that both doublets carry hypercharge y = +1/2, and denote the complex doublet fields by 

<PiW=('?U$), i = l,2. (1) 



For the singlets we assume real fields which we denote by 

<i>i(x), i = l,...,n. (2) 

Of course any complex singlet field may be decomposed into two real singlet fields. 

We employ the bilinear approach for the Higgs doublets (see [16] for details). The Higgs-boson doublets can only 
appear in the form of scalar products flfj with i,j g {1,2} in the potential, ensuring gauge invariance. As shown 
in |16| we may replace the scalar products by the bilinears Kq, K\, K 2 and K3, 

V \i Pl = {K +K s )/2, ip\cp 2 = (K 1 +iK 2 )/2, lf y 2 = {K -K 3 )/2, ^ = (K x - iK 2 )/2. (3) 

The bilinears have to fulfill the conditions 

K > 0, Kl - K\ - K\ - K\ > 0. (4) 

As shown in [TB] the four quantities K a , satisfying Q, parameterize the gauge orbits of the Higgs doublets. The 
advantage of the replacement of the Higgs-boson doublets by the bilinears is that we eliminate the SU(2)l x U(1)y 
gauge degrees of freedom and moreover simplify the potential since the K a are bilinear in the Higgs-boson doublets. 
We thus end up with a potential of the form V(K , K\, K 2 , K3, <f>\, . . . , tf> n ). 



A. Stationary points 

The most general Higgs potential of two Higgs doublets and n real singlets reads 

V(K , Ki,K 2l K 3 , 0i, . . . , <j> n ) =£ a K a + r/ap K a Kp + 

at 4>i + hj cf>i(j)j + c ijk (t>i4>j4>k + dijki 4>i4> 4>k<t>i+ (5) 

&i a (Pi^-a ' Jij a (Pi^Pj^a 

with a,/3 € {0, ...,3} and i,j,k,l £ {l,...,n} and summation convention for repeated indices. Any additional term 
violates SU(2)l X U(1)y invariance respectively renormalizability Note that the coefficients rj a p, bij, c^k, ckjkh an d 
Jij a are symmetric in their greek respectively latin indices. Depending on the number of real singlets n we count 
14 + 5n + 5(' 1 ^ 1 ) + ("J 2 ) + ("J 3 ) coefficients in the most general case. Hence for the case n = 5 we have 219 coefficients 
in the general case. 

In order to determine the stationary points of the Higgs potential d5l) we take the constraint Q into account. 
We distinguish the possible cases of stationary points with respect to the SU(2)^ x U(l)y symmetry breaking be- 
haviour [131 [IB] : 

• Unbroken electroweak gauge symmetry: This is a stationary point with 

K a = K x = K 2 = K 3 = 0. (6) 

A global minimum of this type has vanishing vacuum expectation values for the doublet fields (fl ) . The stationary 



points of this type are found for vanishing Higgs-doublet fields, corresponding to vanishing 
requiring a vanishing gradient with respect to the remaining real singlets: 



ailinears K n and 



V%...,W = 0. (7) 

• Fully broken electroweak gauge symmetry: This is a stationary point with 

K a > 0, K%- K{ - K\ - K\ > 0. (8) 

A global minimum of this type has non- vanishing vacuum-expectation values for the charged components of the 
doublets fields in (Jij, thus gives fully broken SU(2) L x U(1)y ', see [IB]. The stationary points of this type are 
found by requiring a vanishing gradient with respect to all bilinears K a and all singlet fields: 

V V(K , KuK^Ka, &,. ..,&,) = 0. (9) 

In this case the constraints (fSJ) on the bilinears must be checked explicitly for the real solutions. 

• Partially broken electroweak gauge symmetry: This is a stationary point with 

K > 0, Kl - K\ - K\ - K\ = 0. (10) 

For a global minimum of this type follows the desired partial breaking of SU(2)l x U(1)y down to U(l) em ; 
see [IBJ. Using the Lagrange multiplier method, these stationary points are given by the real solutions of the 
system of equations 

V[V{K ,K U K 2 ,K 3 ,4n, . . . ,4>n) - u ■ {Kl - Kl - K\ - K\) ] = 0, 

Kl - Kl - Kl - Kl = 0, 



for (111 



where u is a Lagrange multiplier. The inequality in (10 1 must be checked explicitly for the real solutions found 



If the potential is bounded from below, the global minima will be among the stationary solutions. Plugging the 
solutions into the potential the global minima are the solutions with the lowest potential value. 

The systems of stationary equations, d7|, d9l, (111 are non-linear, multivariate, inhomogeneous systems of polyno- 
mial equations of degree 3. For the case of two doublets and n real singlets we have for the unbroken, fully broken, 
and the partially broken systems sets with n, 4 + n, and 5 + n equations, respectively. The number of indetcrminants 
equals the number of equations. 



III. NUMERICAL POLYNOMIAL HOMOTOPY CONTINUATION METHOD 

Finding the global minima or all stationary points of a given potential is usually a very difficult problem. As outlined 
in the Introduction, a systematic approach is to solve the stationarity systems of equations. While the problem of 
finding all solutions which include minima, maxima and saddle points, compared to only finding minima, seems more 
complicated at first sight, a powerful method called the Grobner basis method can be used to solve the stationary 
equations as long as they are polynomial equations. However, due to the algorithmic complexity issues discussed in 
the Introduction, we are practically restricted to solve only small systems. In particular, going beyond the NMSSM, 
involving two Higgs-boson doublets and one complex singlet (equivalent to two real singlets) , treated in [T3] , appears 
to be difficult to handle with the Grobner basis approach. 

Here, we present a novel numerical method, the numerical polynomial homotopy continuation (NPHC) method, 
which overcomes all the shortcomings of the Grobner basis method. The NPHC method was recently introduced in 
particle theory and statistical mechanics areas in (T7j to find all the so-called Gribov copies of the Landau gauge- 
fixing conditions on the lattice [18 20 and subsequently found many applications in many areas of theoretical physics 
including string phenomenology, lattice field theories, theoretical chemistry, non- linear dynamics, etc. [21H25J . 

The basic strategy behind numerical continuation methods is: first find the solutions of a simple system of equations 
which shares several important features with the given system. Then, in a second step, starting from each of these 
solutions one continues them towards the given system in a systematic way. Homotopy continuation methods have 
been around already for several decades |26[ [27] , but with more recent machinery like the NPHC method used in the 
present article, the method is guaranteed to find all isolated solutions of systems of polynomial equations [TU [TS] . 

Suppose we want to solve a system of m polynomial equations 

P(x)=\ : | i) (.12) 

\Pm(x) 

in the variables x — (x\, . . . , x m ) T , which is known to have only isolated solutions. Then the Bezout's Theorem asserts 
that a system of m polynomial equations in m variables has at most Y\i=i di isolated complex (which obviously include 
real) solutions where di is the degree of the ith polynomial. This upper bound is called the classical Bezout bound, 
and it is known to be sharp for generic systems (i.e., for generic values of the coefficients of the polynomials pi(x)). 
Then, we formulate a homotopy (or, a set of parametric equations) as 

H(x,t) = P(x)(l-t)+jtQ(x), (13) 

where 7 is a complex number and 

Q(x) = : 

\q m (x) 

is again a system of m polynomial equations. Now, varying the parameter t £ [0,1], H can be deformed from the start 
system H(x, 1) = 7 Q(x) at t = 1 into the polynomial system we want to solve, H(x, 0) = P(x) at t — 0. The following 
conditions have to be satisfied in order to guarantee that all solutions of P can be computed from this homotopy: 

1. The solutions of Q(x) = can be computed relatively straightforwardly. 

2. The number of solutions of Q(x) = satisfies the classical Bezout bound for P{x) = as an equality. 

3. The solution set of H (x, t) = for t e (0, 1] consists of a finite number of smooth paths, called homotopy paths, 
which are parameterized by t. 

4. Every isolated solution of H(x,0) — P(x) = 0, including the multiple roots, can be reached by some path 
originating at a solution of H (x, 1) = jQ(x) = 0. 

Satisfying the first two criteria hinges on a suitable choice of the start system Q. Criteria 3 and 4 are guaranteed to 
be satisfied, for generic constants 7 in (l3| [1"4"]. 

The start system Q{x) = can, for example, be taken to be 

(4 1 - 

Q(x)= : j -I). ll.Vi 



= (14) 



where d% is the degree of the i th polynomial of P{x) — 0. The system (15 1 is easy to solve and guarantees that the 
total number of start solutions is Y\i=i di and all solutions are nonsingular. 

Each homotopy path, starting at a solution of Q{x) — at t — 1, is tracked to t = using a path tracking algorithm, 
e.g., Euler predictor and Newton corrector methods. There are a number of freeware packages well-equipped with 
path trackers such as PHCpack [28], HOM4PS2 [25], and Bertini [3D]. We have used Bertini and HOM4PS2 to obtain 
the results in this paper. Tracking the solutions to t = 0, the set of endpoints of these homotopy paths is the set of all 
solutions to P{x) = 0. Since each homotopy path can be tracked independently, NPHC is inherently parallelizable. 
The NPHC method does not suffer from any major computational complexity. It takes floating point coefficients in 
very naturally. 

The set of real solutions can be obtained from the set of complex solutions by considering the imaginary part 
of the solutions (typically, up to a numerical tolerance). We remark that the approach of |31| implemented in 
alphaCertified |32| can be used to certify the reality or non-reality of a nonsingular solution given a numerical 
approximation of the solution. The ability to compute all complex solutions, and thus all real solutions, distinguishes 
the NPHC method from most other methods. 

IV. MINIMIZATION OF THE HIGGS POTENTIAL 
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Figure 1: Higgs potential values for the stationary points with 2 Higgs-boson doublets and 2 Higgs-boson singlets. The 
parameter tan(/3) is varied in the range tan(/3) = 1, . . . , 50. All other parameters are fixed as described in the text. The 
different types of stationary points are shown with respect to the electroweak symmetry breaking behavior. The filled dots 
denote the unbroken, the empty squares the fully broken solutions and the triangles the partially broken solutions, respectively. 



Only the full green triangles correspond to stationary points with vevs as given in (18 1. The stationary point with the lowest 
potential value is the global minimum. 

We now want to apply the NPHC method to the minimization problem of the Higgs potential (pi, i.e., we want 
to solve the stationarity systems of equations d7| , d9l , and (111. We note that the solution with the correct elec- 



troweak symmetry breaking has to come from the set (111. We start with numerically fixing the vacuum-expectation 



values (vevs) v, tan(/3) and v\. with i = 1, ...,n with n being the number of real singlets. As usual tan(/3) = ^2/^1 
denotes the ratio of the vevs of the two Higgs-boson doublets, v 2 — v\ + v\ and v[ are the vevs of the real singlets: 

1 fQ\ , , 1 nl8 f0\ ,,, 1 



w = 7iU • M = T2 e U • <« = 7T' (16) 
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Figure 2: Similar to Fig. [Tithe Higgs potential values for the stationary points with two Higgs-boson doublets but with 3 
Higgs-boson singlets. 

Except for £ a , a = 0, ..., 3 and Oj, i = 1, ..., n, we first fix all coefficients in (pi numerically. The parameters £ Q and a, 
are then fixed by employing the tadpole conditions, 

(Vy( J ftro ) ir i) ^2,^3,01,-.-,^n))=O. (17) 

In this way we gain solutions corresponding to the correct electroweak symmetry breaking behavior with the desired 
vevs. 

For the case considered here with two Higgs-boson doublets and up to five real singlets we choose the vevs 

v = 246 GeV, v[ = 100 GeV, v' 2 = 150 GeV, v' 3 = 200 GeV, v' A = 250 GeV, v' 5 = 300 GeV, (18) 



a vanishing phase 9 of the second doublet and vary the ratio of the doublet vevs in the range tan(/3) = 1, . . . , 50. All 
coefficients in (51, except for £ Q and aj, that is, 7700, 7701, ■ • ■ , /nn3 are fixed in the order they appear in the potential. 
As a numerical example we assign the values 0.1, 0.2, 0.3, 0.4, 0.5, 0.1, . . . ,0.5, 0.1 ... in turn. Note the symmetry 
of the coefficients, that is, only ordered indices appear, for instance 621 is replaced by 612 in the implicit summation 
in ((5| before the values are assigned. 

With all parameters numerically fixed we solve the systems of equations ffi\, ||9J, and (11 1. As discussed in Sec. Ill 



we employ the NPHC method for the solution of the stationary equations. We note that on a regular desktop machine 
we could solve the systems for n = 10 in around 16 hours for a fixed set of parameters. However, since we want to 
scan over a large range of parameter values, due to the available limited computer resources, we restrict ourselves to 
smaller systems, i.e., n — 2, . . . , 5. The largest of these systems can be solved in around 4 hours for a given set of 
parameters on the desktop machine. 

For two Higgs-boson doublets and n — 2, ... ,5 real singlets we show in Table |T] the number of complex stationary 
solutions with respect to the fully, unbroken and partially broken cases. 

From all complex solutions the real solutions are sorted out (practically we require all imaginary parts to be smaller 
than 10 -8 ). The number of real solutions is typically much smaller and depends on the chosen parameters. The fact 
that the number of complex solutions is constant with varied parameters is a good indication that no solutions are 
missed. For the fully broken set subsequently the conditions d8l are checked. Similarly, for the partially broken case, 



the solutions have to fulfill (10) 



The potential values for all stationary solutions passing all checks are plotted in Figs. [TJ [2] [3] [4] for 2, 3, 4, 5 real 
singlets, respectively. The parameter tan(/3) is varied in the range tan(/3) = 1, . . . , 50. Since the potential is bounded 
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Figure 3: Similar to Fig. [Tithe Higgs potential values for the stationary points with two Higgs-boson doublets but with 4 
Higgs-boson singlets. 

n unbroken fully broken partially broken 
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54 
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27 


27 


162 
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75 


81 


486 
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225 


243 


1458 



Table I: Number of complex solutions found for the most general Higgs potential with two doublets and n singlets. The number 
of solutions are given with respect to the systems of equations ¥j\ (unbroken case), (f9j (fully broken case), and ( |11[ ) (partially 
broken case), respectively. 



from below the stationary solutions with the lowest potential value are the global minima. The stationary solutions 
are denoted by circles, squares, and triangles corresponding to unbroken, fully broken, and partially broken solutions, 
respectively. In case of the partially broken solutions empty triangles show solutions which do not give the required 
vevs ( [181 ) > whereas the full triangles show the solutions which yield to the required vevs (with a precision of 8 digits) . 
In order to check whether the partially broken solutions give the desired vevs ( 18 1 we calculate the vevs from the 
solutions via (16), that is, 



^2K =V, yj 






(19) 



We find a very rich structure of stationary points. Obviously the global minimum, that is the vacuum, not always 
correspond to the physically desired solution. As can be seen in Fig. IT] in the case of two real singlets, the global 
minimum breaks clectroweak gauge symmetry partially, however yields undesired vevs for small tan(/3) values. 

Let us note that quite generically we find nearly degenerate potential values for different stationary solutions. 
Note that even when the potential values are nearly degenerate the solutions may be located at very different field 
values. For instance in the case of two Higgs-boson doublets and two singlets we find for tan(/3) = 10 the two deepest 
stationary solutions as shown explicitly in Tab. [IT] In this numerical example only the solution 1 corresponds to the 
desired vevs (18 1, however solution 2 is the global minimum. This behavior of nearly degenerate vacua was already 



mentioned for the case of the NMSSM in [13 which we confirm for more general cases. 
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Figure 4: Similar to Fig. [T] the Higgs potential values for the stationary points with two Higgs-boson doublets but with 5 
Higgs-boson singlets. 



K K x 



K, 



K:>, 



V 



solution f: 30257 5991.6 0.0046390 -29658 70.710 106.06 -1.0334- 10~ 8 -1.5045- 10 9 
solution 2: 34878 8085.9 -4701.0 -33600 55.630 123.54 -0.0036596 -1.5051 - 10 9 



Table II: The two deepest solutions explicitly for the case of two singlets and tan(/3) = 10. Solution 1 corresponds to the desired 
vevs, whereas solution 2 gives wrong vevs. 



V. CONCLUSIONS 

The extended Higgs models have attracted a lot of attention recently, in particular with the on-going experiments 
at the Large Hadron Collider. The main technical difficulty in studying these models is that accurately finding the 
global minimum of the corresponding Higgs potential is highly non-trivial due to its non-linear nature. Whereas the 
Higgs potential of the Standard Model is rather simple, this is in general not longer true for extended potentials. The 
solution of the stationarity equations allows to identify the vacuum, as long as the potential is bounded from below. 
In this paper, we have applied the numerical polynomial homotopy continuation method to solve the non-linear, 
multivariate systems of polynomial stationary equations for rather involved systems of Higgs potentials. In contrast 
to most minimization methods this method guarantees the detection of the global minimum given by the stationary 
solution with the lowest potential value. 

We have shown that the case of the most general Higgs potential with two complex doublets and five real singlets 
is solvable via numerical polynomial homotopy continuation. The results show a rich structure of stationary solutions 
with local maxima, saddle points and minima. In this way, we have shown that a large set of parameter values can be 
excluded from the physically viable regions. Among the solutions we typically find solutions with nearly degenerate 
potential values, however far away in terms of the Higgs-boson fields. Let us note that the numerical polynomial 
homotopy continuation method allows to go gar beyond the limits of two doublets and five singlets. 
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